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Abstract 

Within a general cluster framework, we discuss the loop-algorithm, a new type 
of cluster algorithm that reduces critical slowing down in vertex models and 
in quantum spin systems. We cover the example of the 6- vertex model in 
detail. For the F-model, we present numerical results that demonstrate the 
effectiveness of the loop algorithm. We discuss how to modify the original 
algorithm for some more complicated situations, especially for quantum spin 
systems in one and two dimensions. 



*To appear in Computer Simulations in Condensed Matter Physics VI, ed. D.P. Landau, K.K. Mon, 
and H.B. Schxittler, Springer Verlag, Heidelberg, Berlin, 1993. 



1 Introduction 



For Monte Carlo simulations of many interesting physical situations, critical slowing down 
is a major problem. Standard simulation algorithms employ local update procedures like 
the Metropolis and the heat bath algorithm. With local updates, "information" travels 
slowly, like a random walk. If the relevant length scale is the correlation length ^, the 
number of updates necessary to decorrelate large regions, i.e. the autocorrelation time r, 
grows like 

rcxr, (1) 

where z 2 for local updates, as suggested by the random walk analogy. 

The way out of this problem is to employ nonlocal updates. The challenge is to devise 
algorithms that are nonlocal and still satisfy detailed balance. Multigrid algorithms are one 
possible approach. In this paper we shall focus on cluster algorithms; for a nonexhaustive 
selection of references see [0, ^, ^, ^ |6| . 

The first cluster algorithm was invented by Swendsen and Wang (SW) [j^ for the case 
of the Ising spin model. The basic idea is to perform moves that significantly change the 
Peierls contours characterizing a configuration. As the size of Peierls contours is, typically, 
anything up to the order of the correlation length, critical slowing down may be eliminated 
completely or at least partially by this approach. The SW algorithm has been modified 
and generalized for other spin systems, mostly with two spin interactions |2|, |^, |5[. Notice 
that for these systems clusters are connected regions of spins, with the same dimensionality 
as the underlying lattice. A few generalizations along different lines were also done [Q, ^. 

Recently [0, |8|, we introduced cluster algorithms for vertex models and quantum spin 
systems, which are the first cluster algorithms for models with constraints. While 0| is an 
adaptation of the valleys-to mountains-reflections (VMR) algorithm , originally devised 
for solid-on-solid models, the loop algorithm introduced in ^ does not resemble any 
existing scheme. 

In this paper we discuss the loop algorithm in detail. In vertex models [10] the dynam- 
ical variables are localized on bonds, and the interaction is between all bonds meeting at 
a vertex. Furthermore there are constraints on the possible bond variable values around a 
vertex. Our scheme is devised such as to take into account the constraints automatically, 
and to allow a simple way to construct the clusters. The clusters here are not connected 
regions of spins, but instead closed paths (loops) of bonds. 

In what follows we flrst comment on the SW algorithm for spin systems. Then we 
discuss the general cluster formalism of Kandel and Domany Q , treating the SW algorithm 
as an example. Next we define the 6-vertex model, and, as a special case, the F model. We 
then introduce the loop algorithm, and show how to formulate it for the complete 6-vertex 
model. The optimization of the algorithm is discussed in a separate section. For the special 
case of the F-model, we describe how the algorithm particularizes, and then we present our 
very successful numerical tests of the loop algorithm. We also show how to apply the loop 
algorithm to simulations of one and two dimensional quantum spin systems |jll|], like e.g. 
the xxz quantum spin chain, and the spin ^ Heisenberg antiferromagnet and ferromagnet 
in two dimensions. Finally we comment on further generalizations and applications. 
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2 Some comments on the Swendsen-Wang algorithm 



Here we present a somewhat unusual viewpoint on the SW algorithm, in order to better 
understand what is new in the loop algorithm. Let us look at a spin system, like the Ising 
model, with variables Sx living on sites x, and a nearest neighbor Hamiltonian Hi^Sx^Sy^. 
The partition function of the model to be simulated is Z = exp (— X]<a;j/> H{sx, Sy)), 
where < xy > is a pair of sites. 

Consider update proposals (flips) Sx —>■ s'x, such that H{sx,Sy) = H(s'x,Sy). Then it 
follows that updating at the same time all spins in some "cluster" of spins, we will only 
change the value of the Hamiltonian at the boundary of that cluster, not inside it. 

In order to satisfy detailed balance, we have to choose clusters with an appropriate 
probability. The SW algorithm amounts to making a Metropolis decision for each bond 
< xy >, namely whether to accept the change in H from H{sx,Sy) to H(s'x,Sy). If 
accepted, a bond is called "deleted", otherwise it is called "frozen". Clusters are then sets 
of sites connected by frozen bonds. Note that if deleted, a bond may be at the boundary 
of a cluster, but need not. 

Finally, an update is performed by finding all the clusters in a given configuration and 
then flipping each cluster with 50% probability Q. In Wolff's single cluster variant |^, ^, 
which is dynamically more favourable, we construct one cluster starting from a randomly 
chosen initial site, and then flip it with 100% probability. 

A technical remark: The Swendsen-Wang algorithm can be vectorized and parallelized 
|12| . The difficult task is to identify the clusters, which is the same task as e.g. image 
component labeling. For the single cluster variant, vectorization is the most efficient 
approach |p^ . 



3 The Kandel - Domany framework. 

Cluster algorithms are conveniently described in the general framework of Kandel and 
Domany [^. Let us consider the partition function 

Z = ^exp(-yH), (2) 

u 

where u are the configurations to be summed over. We shall call the function V the 
"interaction". Let us also define a set of new interactions Vi (the index i numbers the 
modified interactions). Assume that during a Monte Carlo simulation we arrived at a 
given configuration u. We choose a new configuration in a two step procedure. The first 
step is to replace V with one of the Vi. For a given i, Vi is chosen with probability Pi{u). 
The Pi{u) satisfy: 

Pi{u) = exp{V{u) -Vi{u) +Ci), 
T,iPi{u) = 1, 

where q are constants independent of u. The second step is to update the configuration 
by employing a procedure that satisfies detailed balance for the chosen Vi. The combined 
procedure satisfies detailed balance for the original interaction V [Q]. 
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In many cases, the interaction y is a sum over local functions H^, where c typically are 
cells of the lattice (like bonds in the spin system case, sites for vertex models, plaquettes for 
gauge theories). More generally, H'^ can contain part (or all) of the interactions in some 
neighborhood of the cell c. We can choose the modified interactions H^ic for each 
independently (i^ numbers the possible new interactions for the cell c). The probabilities 
for this choice obey eq. (^), with V replaced by W^, and i by i'^. The total modified 
interaction Vi will then be the sum of all the H'^ic {i is now a multiindex). When clear 
from the context, we shall drop the index c. 

Take the Ising model as an example. The configuration u is comprised of spins Sx = ±1, 
the cells where the interaction is localized are bonds < xy >, and V{u) = —JJ2<xy> ^xSy 
We choose the bonds as our cells c, so we can perform the first step of the Kandel-Domany 
procedure separately for each bond. The original cell interaction is 

H<^y>{s) = -JsxSy . (4) 

Now, let us define two new bond interactions; the first {i'^^y> = 1) is called "freeze", the 
second (i<^f> = 2) is called "delete": 

' 0, 



^<^'^>freeze(^) = I eL ' ^<^-^>delete(^) = • (5) 

Then from (^) we obtain the Swendsen-Wang probabilities if we choose the constants Cj 
that minimize freezing: 

^'deletc^*) = exp(-J(s^.Sy + 1)) = min(l,exp (-2Jsa;Sy)) , 

(6) 

<xy> / \ , <xy> / \ ^ ' 

^'freeze'^'^'' ~ ^ ~ ^'delete*^^-' ' 
which is just the Swendsen-Wang probability. 

4 The 6- Vertex Model 

The 6- vertex model ||l^ is defined on a square lattice. On each bond there is an Ising-like 
variable that is usually represented as an arrow. For example, arrow up or right means 
+1, arrow down or left means —1. At each vertex we have the constraint that there are 
two incoming and two outgoing arrows. In fig. 1 we show the six possible configurations 



at a vertex, numbered as in [10|. The statistical weight of a configuration is given by the 
product over all vertices of the vertex weights p{u). Thus a priori, for each vertex there are 
6 possible weights p{u), u = 1, ...,6. We take the weights to be symmetric under reversal 
of all arrows. Thus, in standard notation [l^], we have: 

p{l)=p{2) = a, 

pi3) = p{4) = 6 , (7) 
p(5) = p(6) = c . 

The 6-vertex model has two types of phase transitions: of Kosterlitz-Thouless type and 
of KDP type |jlO|. A submodel exhibiting the former is the F model, defined by c = 1, 
a = b = exp (—K), K > 0. For the latter transition an example is the KDP model itself, 
defined by a = 1, 6 = c = exp {—K), K > 0. 
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1 2 3 4 5 6 

Figure 1: The six vertex configurations, u = 1,...,6 (using the standard conven- 



tions of 1^10 
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ll-ur ul-lr straight 

Figure 2: The three break-ups of a vertex: ll-ur (lower- left -upper-right), ul-lr 
(upper- left -lower-right), and straight. 



5 The Loop Algorithm 

If we regard the arrows on bonds as a vector field, the constraint at the vertices is a zero- 
divergence condition. Therefore every configuration change can be obtained as a sequence 
of loop-flips. By "loop" we denote an oriented, closed, non-branching (but possibly self- 
intersecting) path of bonds, such that all arrows along the path point in the direction of 
the path. A loop-flip reverses the direction of all arrows along the loop. 

Our cluster algorithm performs precisely such operations, with appropriate probabili- 
ties. It constructs closed paths consisting of one or several loops without common bonds. 
All loops in this path are flipped together. 

We shall construct the path iteratively, following the direction of the arrows. Let the 
bond h be the latest addition to the path. The arrow on h points to a new vertex v. There 
are two outgoing arrows at v, and what we need is a unique prescription for continuing 
the path through v. This is provided by a break-up of the vertex v. In addition to the 
break-up, we have to allow for freezing of v. By choosing suitable probabilities for break-up 
and freezing we shall satisfy detailed balance. 

The break-up operation is defined by splitting v into two pieces, as shown in fig. 2. 
The two pieces are either two corners or two straight lines. On each piece, one of the 
arrows points towards v, while the other one points away from v. Thus we will not allow, 
e.g., the ul-lr break-up for a vertex in the configuration 3. If we break up v, the possible 
new configurations are obtained by flipping (i.e., reversing both arrows of) the two pieces 
independently. On the other hand, if we freeze v, the only possible configuration change 
is to flip all four arrows. 

The break-up and freeze probabilities are conveniently described within the Kandel 
Domany framework. It is sufficient to give them for one vertex, which is in the current 
configuration u. We define 6 new interactions (weight functions) pi, i = 1,...,6, corre- 
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i 


action 




Pi{u) 


1 


freeze 1,2 


1, p{u) = a 
0, else 


qi/a, p[u) = 
a 

(l £il CCl 

dot; 


2 


freeze 3,4 


1, p{u) = h 
0, else 


q2/b, p(u) = 
b 

U, else 


3 


freeze 5,6 


1, p{u) = c 
0, else 


gs/c, p[u) = 
c 

0, else 


4 


11-ur 


0, p{u) = a 

1, else 


0, p(u) = a 

Q4/p{u), 

else 


5 


ul-lr 


0, p{u) = b 

1, else 


0, p{u) = b 
else 


6 


straight 


0, p{u) = c 

1, else 


0, p[u) = c 
else 



Table 1: The new interactions Pi{u) and the probabilities Pi{u) to choose them at a 
vertex in current configuration u. See eq. (0). 



sponding to specific break-up and freeze operations (the labelling of the new interactions 
is completely arbitrary, and the fact that we have six of them is just a coincidence). For 
each vertex in configuration u, we replace with probability Pi{u) the original interaction p 
by the new interaction pi. Equation ^ , i.e. detailed balance and the proper normalization 
of probabilities, require that for every u 

Pi{u) = qi^^4-l , y'pi('u) = 1 , (8) 

where Qi = exp (— q) > are parameters. 

As discussed in Q (see also table 1), freezing is described by introducing one new 
interaction for each different value of p{u). For example, to freeze the value a, we choose 
the interaction pi to be pi{u) = 1 if p{u) = a, and pi{u) = otherwise. In other words, 
when pi is chosen, transitions between u = 1 and u = 2 cost nothing, whereas the vertex 
configurations 3, 4, 5, and 6 are then not allowed. Notice that we denote by u the current 
configuration, and by u the argument of the function pj. 

Each break-up is also described by one new interaction. As an example take the ul-lr 
break-up. It is given by the new interaction number five, with P5(n) = 1 if p{u) = a 
or c, and P5{u) = if p{u) = b. In other words, with the new interaction p^, transitions 
between 1, 2, 5 and 6 cost nothing, while the vertex configurations 3 and 4 are not allowed. 
This corresponds precisely to allowing independent corner flips in a ul-lr break-up (see 
figs. 1,2). 
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Table 1 gives the full list of new weights Pi{u) and probabilities Pi{u) to choose them 
in current configuration u. From (^) we also obtain: 

qi + q5 + q6 = a , 

q2 + q4: + q6 = b , (9) 

93 + 94 + 95 = C . 

Assume now that we have broken or frozen all vertices. Starting from a bond bo, we 
proceed to construct a closed path by moving in the arrow direction. As we move from 
vertex to vertex, we always have a unique way to continue the path. At broken vertices the 
path enters the vertex through one bond and leaves it through another. If the last bond b 
added to the cluster points to a frozen vertex v, the path bifurcates in the directions of the 
two outgoing arrows of v. One of these directions can be considered as belonging to the 
loop we came from, the other one as belonging to a new loop. Since we also have to flip 
the second incoming arrow of v, we are assured that this new loop also closes. The two 
loops have to be flipped together. In general, the zero-divergence condition guarantees 
that all loops will eventually close. 

We have now finished describing the procedure for constructing clusters. In order to 
specify the algorithm completely, we must choose values for the constants qi, and decide 
how the clusters are flipped. The former problem is of utmost importance, and it is the 
object of the next chapter. For the cluster flips, we may use both the Swendsen-Wang 
procedure and the Wolff single cluster flip Q. We choose the latter, i.e. we "grow" a single 
path from a random starting bond bo, and flip it. The break-or- freeze decision is only 
needed for the vertices along the path, so the computer time for one path is proportional 
to the length of that path. 

There are some distinct differences between our loop-clusters and more conventional 
spin-clusters. For spin-clusters, the elementary objects that can be flipped are spins; 
freezing binds them together into clusters. Our closed loops on the other hand may be 
viewed as a part of the boundary of spin-clusters (notice that the boundary of spin clusters 
may contain loops inside loops). It is reasonable to expect that in typical cases, building 
a loop-cluster will cost less work than for a spin-cluster. This is an intrinsic advantage of 
the loop algorithm. 

6 Optimization of free parameters 

We have seen that freezing forces loops to be flipped together. Previous experience with 
cluster algorithms suggests that it is advantageous to be able to flip loops independently, 
as far as possible. We therefore introduce the principle of minimal freezing as a guide 
for choosing the constants g,: we shall minimize the freezing probabilities, given the con- 
straints @ and Qi > 0. In the next chapter we will show that for the case of the F model, 
optimization by minimal freezing does indeed minimize critical slowing down. Here we 
discuss optimization for the 4 phases of the 6-vertex model, usually denoted by capital 
roman numerals [|lO| ]. 
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Let us first look at phase IV, where c > a + b. To minimize the freezing of weight c, 
we have to minimize ^3. From (^), (73 = c — a — 6 + qi + ^2 + 2^6- With Qi > this imphes 
% min = c — a — b. The minimal value of can only be chosen if at the same time we 
set qi = q2 = 0, i.e. minimize (in this case do not allow for) the freezing of the smaller 
weights a and b. The optimized parameters for phase IV are then: 

qi = 0, 92 = 0, qs = c - a - b, 
qA = b, q5 = a, qQ = 0. 

In phase I the situation is technically similar. Here a > b + c, and we minimize freezing 
with qi = a — b — c and q2 = Qs = 0. The same holds for phase II, 5 > a + c, where we 
obtain minimal freezing for q2 = b — a — c and qi = q^ = 0. 

Phase III (the massless phase) is characterized by a, 6, c < ^{a + b + c). Here we can 
set all freezing probabilities to zero. Thus, 

qi = 0, 2^4 = b + c — a , 

92 = 0, 2q5 = c + a-b, (11) 

93 = 0, 2qQ = a + b- c . 



7 Case of the F model 

The F model is obtained from (10) and (11) as the special case a = b = exp {—K) < 1, 
c = 1. It has a Kosterlitz-Thouless transition at Kc = In 2, with a massless phase for 
K<K,. 

How do we choose the parameters qi here ? Symmetry a = b implies qi = 92 and 
94 = 95. We can eliminate freezing of vertices 1, 2, 3, 4, by setting 91 = 92 = 0. In (|9|) this 
leaves one parameter, 93: 

2Q4 = 1-93, .^2) 
2qG = exp (-i^) + 93 - 1 • ^ ^ 

In the massless phase, we can minimize freezing by also setting 93 = 0. In the massive 
phase, qe > limits 93. Thus 

_/ l-2exp(-i^), K>K„ 
93,min - I 0, K<K,. 

Notice that since a = b in the F model, the straight break-up, the freezing of a, and that 
of b are operationally the same thing. 

If we choose 93 = 93 jj^jj^, then for K < Kc vertices of type 5 and 6 are never frozen, 
which has as a consequence that every path consists of a single loop. This loop may 
intersect itself, like in the drawing of the figure "8". For K > Kc on the other hand, 
vertices of type 1, 2, 3 and 4 are never frozen, so we do not continue a path along a 
straight line through any vertex. As — (X) (temperature goes to zero), most vertices 
are of type 5 or 6, and they are almost always frozen. Thus the algorithm basically flips 
between the two degenerate ground states. 
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For the F model we also have a spin-cluster algorithm - the VMR algorithm 0. At 
K = Kc and for 53 = ^3 j^j^, we have a situation where the loop-clusters form parts of the 
boundary of VMR spin-clusters. Since flipping a loop-cluster is not the same as flipping a 
VMR cluster, we expect the two algorithms to have different performances. We found (see 
and the next section) that in units of clusters, the VMR algorithm is more efficient, 
but in work units, which are basically units of CPU time, the loop algorithm wins. At 
where the loop-clusters are not at all related to the boundary of VMR clusters, we 
found the loop algorithm to be more efficient both in units of clusters and in work units, 
with a larger advantage in the latter. 

8 F-model: Performance of the loop algorithm 

We tested the loop algorithm on LxL square lattices with periodic boundary conditions at 
two values of K: at the transition point and at ^Kc, which is deep inside the massless 
phase. We carefully analyzed autocorrelation functions and determined the exponential 
autocorrelation time r. At infinite correlation length, critical slowing down is quantified 
by the relation ([l|) , r oc . 

Local algorithms are slow, with z ~ 2. For comparison, we performed runs with a local 
algorithm that flips arrows around elementary plaquettes with Metropolis probability, and 
indeed found z = 2.2(2) at = Kc- 

In order to make sure that we do observe the slowest mode of the Markov matrix, we 
measured a range of quantities and checked that they exhibit the same r. As in 0, it 
turned out that one can use quantities defined on a sublattice in order to couple strongly to 
the slowest mode. Specifically, we wrote the energy as a sum over two sublattice quantities. 
We shall present more details of this phenomenon elsewhere. Let us however note here 
that for the total energy, the true value of r was not visible within our precision except for 
a weak hint of a long tail in the autocorrelations on the largest lattices we considered. Note 
that as a consequence of this situation, the so-called "integrated autocorrelation time" 
is much smaller than r, and it would be completely misleading to evaluate the algorithm 
based only on its values. 

We shall quote autocorrelation times r in units of "sweeps" defined such that on 
the average each bond is updated once during a sweep. Thus, if r^^ is the autocorrelation 
time in units of clusters, then r = r'^' x (cluster size)/(2L^). Each of our runs consisted 

cl 

of between 50000 and 200000 sweeps. Let us also define z'^^ by r'^' oc , and a cluster 
size exponent c by (cluster size) oc L^. We then have: 

z = - (2 - c) . (14) 

Table 2 shows the autocorrelation time r for the optimal choice = At 
K = ^Kc, deep inside the massless phase, critical slowing down is almost completely 
absent. A fit according to eq. [l| gives z = 0.19(2). The data are also consistent with z = 
and only logarithmic growth. For the cluster size exponent c we obtained c = 1.446(2), 
which points to the clusters being quite fractal (notice that z'^^ = 0.74(2)). At the phase 
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L 


K = Kc 




8 


1.8(1) 


4.9(4) 


16 


3.0(2) 


5.6(2) 


32 


4.9(4) 


6.2(3) 


64 


7.2(7) 


7.4(3) 


128 


15.5(1.5) 


8.3(2) 


256 


20.5(2.0) 




z 


0.71(5) 


0.19(2) 



Table 2: Exponential autocorrelation time r at = (?3,minj and the resulting dynamical 
critical exponent z. 



K 




z 


\Kc 





0.19(2) 




0.10 


1.90(5) 


iKc 


0.20 


> 2.6(4) 


Kc 





0.71(5) 


Kc 


0.05 


0.77(6) 


Kc 


0.10 


0.99(6) 


Kc 


0.20 


> 2.2(1) 



Table 3: Dependence of the dynamical critical exponent z on the parameter q^. We 
use ">" where for our lattice sizes r increases faster than a power of L. 

transition K = Kc we obtained z = 0.71(5), which is still small. The clusters seem to be 
less fractal: c = 1.060(2), so that z""^ = 1.65(5). 

We noted above that at K = Kc and for the optimal choice of gs, the loop-clusters are 
related to the VMR spin-clusters. In 0] we obtained for the VMR algorithm at K = Kc 
the result z"^ = 1.22(2), but we had c = 1.985(4), which left us with z = 1.20(2). In 
this case, although in units of clusters the VMR algorithm is more efficient, the smaller 
dimensionality of the loop-clusters more than make up for this, and in CPU time the loop 
algorithm is more efficient. 

As mentioned, no critical slowing down is visible for the integrated autocorrelation 
time Tint{E) of the total energy. At K = Kc, Ti^t{E) is only 0.80(2) on the largest lattice, 
and we find the dynamical exponent Zint{E) ~ 0.20(2). At K = ^Kc, Tint{E) is 1-1(1) on 
all lattice sizes, so Zuit{E) is zero. 

What happens for non-minimal values of ? Table 3 shows our results on the depen- 
dence of z on z rapidly increases as moves away from g3niin- This effect seems to 
be stronger at ^Kc than at Kc- We thus see that the optimal value of q^ indeed produces 
the best results, as conjectured from our principle of least possible freezing. 
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In the massive phase close to Kc^ we expect that z{Kc) will determine the behaviour 
of r in a similar way as in ref. Q. To confirm this, a finite size scaling analysis of r is 
required. In order to complete the study of the loop algorithm's performance, we should 
also investigate it at the KDV transition. 

In summary, the loop algorithm strongly reduces critical slowing down, from z = 2.2(2) 
for the Metropolis algorithm, down to z = 0.71(5) at and z = 0.19(2) at i^Kc deep 
inside the massless phase. For the integrated autocorrelation time of the total energy, no 
critical slowing down is visible in either case. 
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9 Quantum Spin Systems 



Particularly promising is the possibility of accelerating Quantum Monte Carlo simulations, 
since quantum spin systems in one and two dimensions can be mapped into vertex models 
in 1 + 1 and 2 + 1 dimensions via the Trotter formula and suitable splittings of the 
Hamiltonian |jll[. The simplest example is the spin ^ xxz quantum chain, which is mapped 
directly into the 6-vertex model. For higher spins, more complicated vertex models result 
(e.g. 19- vertex model for spin one). 

For (2 + 1) dimensions, different splittings of the Hamiltonian lead to quite different 
vertex models, in particular on quite different lattice types |11]. For example, in the case 
of spin ^ we can choose between 6-vertex models on a quite complicated 2 + 1 dimensional 
lattice, and models on a bcc lattice, with 8 bonds and a large number of configurations 
per vertex. 

For the simulation of the 2- dimensional Heisenberg antiferromagnet and ferromagnet 
using the former splitting, all relevant formulas have been worked out in the present 
paper. Actually, the low temperature properties of the antiferromagnet have recently 
been investigated by Wiese and Ying [0] using our algorithm. Their calculation is, in 
our opinion, the first high quality verification of the magnon picture for the low lying 
excitations. In particular, this excludes to a much higher degree of confidence than before 
the speculation (some years ago widespread) that the model had a nonzero mass gap. 

Notice that, similarly to other cluster algorithms j^, it is straightforward to define 
improved observables. The investigation |14| in fact uses them. 

Let us also remark that the loop algorithm can easily change global properties like 
the number of world lines or the winding number (see |11]). Thus it is well suited for 
simulations in the grand canonical ensemble. Last, but not least, the loop algorithm also 
opens up a new avenue for taming the notorious fermion sign problem [|l5| . 



10 Conclusions 

We have presented a new type of cluster algorithm. It flips closed paths of bonds in vertex 
models. Constraints are automatically satisfled. We have successfully tested our algorithm 
for the F model and found remarkably small dynamical critical exponents. 

There are many promising and straightforward applications of our approach, to other 
vertex models, and to 1+1 and 2+1 dimensional quantum spin systems. Investigations 
of such systems are in progress. In particular, we believe that our generalizations of the 
freeze-delete scheme can be adapted for other models like the 8-vertex model. 

Already, the loop algorithm has found important applications in the study of the 2- 
dimensional Heisenberg antiferromagnet. 
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